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A  model  is  presented  that  describes  the  main  physical  phenomena  affecting  in  the  performance  of  a  Solid- 
Oxide  Fuel  Cell  (SOFC).  The  implementation  of  the  model  uses  an  in-house  algorithm  in  a  computational 
fluid-dynamics  (CFD)  framework  that  may  be  used  to  optimize  the  SOFC  operational  parameters.  The 
physical  phenomena  considered  in  the  model  are:  (i)  mass  conservation:  multicomponent  and  multi¬ 
modal  mass  transfer  in  gas  channels  and  electrodes  (convection,  ordinary  diffusion,  Knudsen  diffusion); 
(ii)  momentum  conservation  in  the  gas  channels  and  electrodes;  (iii)  energy  conservation:  coupled  heat 
transfer  across  the  whole  cell  (gas  channels,  electrodes  and  electrolyte);  (iv)  electrochemistry:  half¬ 
reactions  are  considered  to  take  place  at  the  electrode-electrolyte  interfaces,  and  activation  losses  are 
computed  using  the  general  version  of  the  Butler-Volmer  equation.  The  main  features  of  this  CFD  tool 
are:  (i)  it  allows  the  prediction  of  the  characteristic  (. I-V )  curve  of  an  H2-fed  cell;  (ii)  it  is  suitable  for 
both  tubular  and  planar  cells;  (iii)  it  has  been  implemented  using  OpenFOAM-1.5-dev,  an  open-source 
CFD-platform  based  on  the  Finite  Volume  Method.The  numerical  results  are  validated  with  published 
experimental  I-V  curves  for  a  hydrogen-fed  anode-supported  micro-tubular  SOFC,  and  a  numerical  anal¬ 
ysis  of  the  influence  of  different  operation  conditions  on  the  temperature  distribution  is  performed  to 
procure  a  better  understanding  of  the  heat  management  of  the  cell. 

©  201 1  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  environmental  concerns  and  the  geopolitical  consequences 
of  the  use  of  fossil  fuels  have  prompted,  in  the  last  decades, 
the  need  for  new  and  cleaner  energy  technologies.  Hence,  solid 
oxide  fuel  cells  (SOFC)  are  increasingly  regarded  as  a  future 
environmental-friendly  power-generation  technology.  Although 
some  pre-commercial  prototypes  are  already  available  [1,2],  a 
wider  use  of  SOFC  still  requires  a  research  effort  to  reduce  the 
high  costs,  to  increase  the  durability  and  to  improve  their  start¬ 
up  performance.  To  reach  these  goals,  current  research  is  focused 
on  the  study  of  new  cell  materials  and  structures,  which  enable 
good  cell  performance  at  lower  operating  temperatures.  The  anode- 
supported  microtubular  SOFC  is  one  of  those  promising  new  cell 
structures  [3,5,6],  since  it  offers  high  thermal  shock  resistance, 
rapid  startability,  lower  operating  temperatures,  higher  power 
densities  and  simpler  seal  requirements.  However,  mass  transport 
and  heat  management  become  critical  during  the  operation  of  such 
micro-tubular  cells,  since  the  high  power  density  enhances  Joule 
heating  and  the  thick  anode  hinders  the  supply  of  the  reactants 
to  the  reaction  sites.  Since  the  species  and  temperature  distribu¬ 
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tions  within  the  cell  are  not  easy  to  measure,  numerical  simulations 
are  being  increasingly  employed  to  understand  the  experimental 
evidence  and  to  steer  the  cell  optimization  [4]. 

The  modeling  by  the  authors  of  mass-transfer  in  SOFCs  was 
already  addressed  in  a  previous  paper  [7,8],  where  multimodal 
mass-transfer  was  thoroughly  studied  under  isothermal  condi¬ 
tions.  In  this  paper,  heat  transfer  issues  in  SOFCs  are  explored.  It 
will  be  shown  that  the  inlet  temperatures  of  the  gases  play  a  crit¬ 
ical  role  in  the  cell  performance,  in  particular  when  the  feeding 
velocities  are  large  and  the  gases  are  not  preheated.  In  these  condi¬ 
tions,  the  convective  cooling  of  the  cell  may  lead  to  a  degradation 
of  performance  that  may  be  wrongly  attributed. 

To  conduct  the  analysis,  the  pre-existing  model  is  extended  with 
heat-transfer  capabilities.  Hence,  the  model  accounts  not  only  for 
the  mass  and  momentum  conservation  in  the  gas  channels  and 
electrodes  and  the  electrochemistry;  but  also  for  the  coupled  heat 
transfer  across  the  whole  cell.  Recent  reviews  on  SOFC  model¬ 
ing  [9,10],  indicate  that  models  with  heat-transfer  features  have 
indeed  been  reported  in  the  past;  among  these,  the  authors  would 
like  to  highlight  the  comprehensive  modeling  presented  by  Ser- 
incan  et  al.  [27].  The  novelty  introduced  by  the  present  work  is 
multifacted.  Mathematically,  a  new  formulation  to  solve  the  cou¬ 
pled  heat  transfer  within  the  cell  is  presented  (see  Section  2), 
where  the  sensible-enthalpy  conservation-equation  is  derived  in 
terms  of  temperature  avoiding  the  use  of  constant  thermodynamic 
properties.  The  temperature-dependence  of  the  transport  and  ther- 
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Nomenclature 

a  Constant,  Eq.  (37) 

A  Area(m2) 

b  Constant,  Eq.  (37) 

B0  Permeability  of  the  porous-medium  (m2 ) 

Cp  Specific  heat  at  constant  pressure  of  the  fluid 

(m2  s-2  K-1 ) 

Dap  Binary  diffusion  coefficient  of  species  a  in  species  /3 

(m2  s-1) 

Dam  Diffusion  coefficient  of  species  a  in  the  gas  mixture 

(m2  s-1) 

D®  p  Effective  binary  diffusion  coefficient  of  species  a  in 

species/3  (m2  s-1) 

Effective  Knudsen  diffusion  coefficient  of  species  a 

(m2  s-1) 

E  Energy  (kg  m2  s-2  kmol-1 ) 

E°  Standard  electrochemical  cell  voltage  (V) 

Eb  Blackbody  emissive  power 

/  Body  forces  (ms-2) 

F  Faraday’s  constant  (A  s  kmol-1 ) 

Fi_j  View  factor  between  surface  i  and  surface  j 
h  Sensible  enthalpy  of  the  fluid  (m2  s-2) 

ha  Sensible  enthalpy  of  the  species  a  in  the  fluid 
(m2  s-2) 

H0  Incident  radiation  entering  or  leaving  the  enclosure 
through  an  opening  (kgs-3) 
i  Unit  vector  in  the  direction  of  x-axis 

I  Current  density  (A  m-2 ) 

J0  Exchange  current  density  (A  m-2 ) 

ja  Mass  diffusive  flux  of  species  a  (kg  m-2  s-1 ) 

k b  Boltzmann’s  constant  (kg  m2  s-2 1<-1 ) 
l  Thickness  (m) 

n  Number  of  electrons  to  convert  a  single  molecule  of 

species  a 

h  Normal  unit-vector 

Na  Total  molar  flux  of  species  a  (kmol  m-2  s-1 ) 
p  Pressure  (kgm-1  s-2) 

Pa  Partial  pressure  of  species  a  (kgm-1  s-2) 

P  Power  (kgm2  s-3) 

q  Energy  flux  (kg  s-3 ) 

Q  Volumetric  heat  sources,  Joule  heat  (kg  m-1  s-3 ) 

Qin  Volumetric  inlet  flow  (m3  s-1 ) 

r  Vector  from  channel  surface  to  enclosure  (m) 

R  Ideal  gas  constant  (kg  m2  s-2  kmol-1  K-1 ) 

S  Volumetric  mass  source  (kg  m-3  s-1 ) 

Sa  Volumetric  mass  source  of  species  a  (kg  m-3  s-1 ) 

Sx  Sutherland-law  parameter  (K) 

5^  Sutherland-law  parameter  (K) 

T  Temperature  (I<) 

Tox  Sutherland-law  parameter  (K) 

Tqjjl  Sutherland-law  parameter  (K) 

v  Fluid  velocity  (ms-1) 

V  Voltage  (V) 

Wa  Molecular  weight  of  species  a  (kg  kmol-1 ) 
xa  Molar  fraction  of  species  a  in  the  gas  mixture 
ya  Mass  fraction  of  the  species  a 

Greek  symbols 

a  Gas  mixture  species 

a  Backward  transfer  coefficient 

a  Forward  transfer  coefficient 

/3  Gas  mixture  species 

y  Pre-exponential  coefficient  (A  m-2 ) 


r*  Dusty-gas  model  parameter  (kg-1  s  kmol) 

s  Porosity 

sa  Characteristic  Lennard-Jones  energies  of  species  a 

(kg  m2  s-2) 

£rad  Emissivity  of  the  surface 

rj  Overpotential  (V) 

9  Angle 

A  Thermal  conductivity  of  the  fluid  (kg  m  s-3 1<-1 ) 

As  Thermal  conductivity  of  the  solid  (kg  m  s-3  I<-1 ) 

Aeff  Effective  thermal  conductivity  of  the  fluid  and 
porous  medium  as  a  continuum  (kg  m  s-3  K-1 ) 

A0  Sutherland-law  parameter  (kg  m  s-3  K-1 ) 
pi  Viscosity  of  the  fluid  (kgm-1  s-1) 

fia  Viscosity  of  species  a  (kg  m-1  s-1 ) 

/jlo  Sutherland-law  parameter  (kg  m-1  s-1 ) 

p  Fluid  density  (kg  m3 ) 

a  Anionic  conductivity  (A  V-1  m-1 ) 

crap  Collision  diameter  (Lennard-Jones  12-6  potential 

model)  (A) 

<ja  Characteristic  length  of  species  a  (A) 

ctS-b  Stefan-Boltzmann  constant  (kg  s-3 1<-4) 

r  Tortuosity  factor 

f  Viscous  stress  tensor  (kgm-1  s-2) 

u£*  Dusty-gas  model  parameter  (kg-1  s  kmol) 

uJJJ*  Dusty-gas  model  parameter  (kg-1  s  kmol) 

</>v  Viscous  dissipation  (kgm-1  s-3) 

zu  Heat  of  reaction  (kg  m-1  s-3 ) 

zu'  Molar  reaction  heat  release  (kg  m2  s-2  kmol-1 ) 

£2d  Collision  integral 

Subscripts 
a  Anode 

act  Activation 

c  Cathode 

Cl  Channel  interface  coupled  to  the  electrode 

con  Concentration 

e  Electrolyte 

El  Electrode  interface  coupled  to  the  channel 

fuel  Fuel 

i  Index  for  channel  isothermal  differential  surfaces 

in  Channel  inlet 

j  Index  for  enclosure  isothermal  differential  surfaces 

new  New  iteration 

ocv  Open  circuit  voltage 

ohm  Ohmic 

out  Channel  outlet 

oven  Heated  furnace  wall 

rad  Radiation 

ref  Reference 

rw  Reaction  wall,  active  electrode-electrolyte  interface 


modynamic  properties  of  the  gases  and  solid  components  of  the  cell 
is  fully  considered,  enhancing  the  accuracy  of  the  temperature  cal¬ 
culation.  From  the  point  of  view  of  Computational  Fluid  Dynamics 
(CFD),  unlike  most  of  the  previous  SOFC  models,  which  are  based  on 
ad-hoc  extensions  to  commercial  codes,  the  numerical  algorithm 
presented  in  this  paper  (Section  3)  has  been  fully  implemented  by 
the  authors  in  OpenFOAM,  an  open-source  CFD  platform.  Finally, 
in  Section  4,  the  model  has  been  exploited  from  a  scientific  point 
of  view  to  get  a  better  understanding  of  the  heat  management  in 
SOFCs.  In  Section  4,  the  validity  of  the  model  is  proved  by  com¬ 
parison  of  its  results  against  experimental  data  [5];  and  it  is  then 
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used  to  study  the  cooling  effect  of  the  feeding  fuel  on  the  cell  Finally,  energy  conservation  in  the  system  is  considered  through 

performance.  the  sensible-enthalpy  equation: 


2.  Mathematical  model 

In  this  section,  the  mathematical  model  that  describes  the 
steady-state  operation  of  a  solid  oxide  fuel  cell  is  described.  This 
model  is  based  on  a  previous  one  by  the  same  authors  for  mass- 
transfer  in  SOFCs  [7,8],  which  has  been  extended  to  consider  other 
relevant  phenomena  involved  in  the  SOFC  operation:  heat  trans¬ 
fer  and  electrochemistry.  Due  to  multiphysics  nature  of  SOFCs,  the 
model  consists  of  four  subsets  of  equations:  (i)  channel  model;  (ii) 
electrode  model;  (iii)  electrolyte  model;  and  (iv)  electrochemistry 
model. 

2.1.  Channel  model 

The  set  of  equations  solved  to  model  both  channels  are  continu¬ 
ity,  momentum,  species  and  enthalpy  conservation-equations. 

The  continuity  equation  is: 

V-(p5)  =  S  (1) 

where  p  is  the  fluid  density,  v  is  the  fluid  velocity  in  the  channel  and 
S  is  the  volumetric  mass  source  term.  The  momentum-conservation 
equations  are: 

V  •  {pvv)  -  V  •  x  —  -Vp  (2) 

where  r  is  the  viscous  stress  tensor  and  p  is  the  pressure.  The 
equation  for  the  conservation  of  the  chemical  species  a  is  written 
as: 

v  •  (pyav)  +  v  •  ja  =  Sa  (3) 


V  •  (phv)  +  V  •  (3)  =  V  ■  Vp  +  0v  +  Q  +  P^JafaVg  +  7JJ  (8) 

a 

where  h  is  the  sensible  enthalpy  of  the  fluid  (h  =  ^ayaha),  q  is 
the  energy  flux,  is  the  viscous  dissipation,  Q.  represents  the  vol¬ 
umetric  heat  sources,  fa  accounts  for  the  body  forces  and  m  for 
the  heat  of  reaction.  The  following  assumptions  apply  to  the  study 
of  the  SOFC  operation:  the  effect  of  body  forces  has  no  relevance 
on  the  system  energy  {p^2ayoJaVa  ~  0);  energy  dissipation  due  to 
viscous  forces  in  gas  flows  under  laminar  regime  may  be  neglected 
(0V  ~  0);  energy  sources  due  to  compressibility  effects  are  ignored, 
since  large  pressure  differences  are  not  expected  inside  the  cell 
(v-  Vp  «  0);  and  the  heat  flux  accounts  for  the  heat  conduction 
and  the  heat  flux  due  to  species  diffusion  with  different  enthalpies 
( q  =  -AVI  +  'Ejaha).  Then,  Eq.  (8)  is  simplified  to: 

V  ■  (phv)  -  V  •  (XVT)  +  V  ■  =Q-  +  m  (9) 

In  the  SOFC  channels,  the  reaction  heat  release  is  zero  (zzr  =  0), 
because  the  reaction  sites  are  located  in  the  electrodes.  According 
to  previous  studies  [20,21  ],  the  gases  in  the  SOFC  can  be  regarded  as 
non-participating  media  for  radiation.  In  the  cathodic  channel,  the 
air  consists  of  simple  non-polar  molecules,  transparent  to  radiation. 
In  the  anodic  channel,  the  radiation  of  the  participating  species  has 
been  shown  to  have  a  negligible  effect  on  the  SOFC  performance, 
under  normal  SOFC  operation  conditions.  Thus  a  surface-to-surface 
radiation  model  has  been  used  and  will  be  described  later;  but,  for 
the  fluids  in  the  channels,  there  are  no  volumetric  heat  sources 
due  to  radiation  and  Q  =  0  above.  The  final-heat  transfer  equation  is 
therefore: 


where  ya  is  the  mass  fraction  of  the  species  aja  is  the  mass  diffusive 
flux  of  species  a ,  and  Sa  is  the  volumetric  source  term  for  the  species 
a.  The  multicomponent  diffusion  is  modelled  using  a  consistent 
effective  binary  diffusion  method  [11]: 

ja  = -pDamVya+yap^2DpmVyp  (4) 

VyS 


where  Dam  is  the  diffusion  coefficient  of  species  a  in  the  gas  mix¬ 
ture,  given  by: 


D 


am 


1  —  Xa 


(5) 


Here,  xa  is  the  molar  fraction  of  species  a  and  Dap  is  the  binary 
diffusion  coefficient  of  species  a  in  species  (B,  modelled  as  [13]: 


Dap  —  2.628  x  10“3 


y/T3(Wa  +  W^2WaW^) 

CTaP^DP 


(6) 


where  Wa  and  Wp  are  the  molecular  weights  of  species  a  and 
p  respectively;  crap  is  the  collision  diameter  (Lennard-Jones  12-6 
potential  model)  given  by  crap  =  (cra  +  crp)/2,  with  cra  and  crp  being 
the  characteristic  length  of  species  a  and  |3;  and  is  the  collision 
integral: 

1.06036  0.1930  1.03587 

D  -  (T^o.  15610  +  exp(0.47635r*)  +  exp(  1.529967*) 

1.76474 

+  exp(3.89411T*)  (  ’ 


Flere,  T*  =  /<BT(£a£p)_a5,  where  /<B  is  the  Boltzmann’s  constant 
and  £a  and  £p  are  the  characteristic  Lennard-Jones  energies  of 
species  a  and  p  respectively  [14].  Data  for  £a,  £p,  <raand  crp  are 
taken  from  [15]. 


V  •  (phv)  -  V  •  (AVI)  +  V  •  \^2jJaha)j  =  0  (10) 

Considering  the  definition  of  the  sensible  enthalpy  (dh  =  CpdT),  Eq. 
(10)  may  be  rewritten  as: 

V-(pCpvT)-V-(XVT)  =  TV-(Cppv)-V-  (j^Oa/ia))  (H) 

where  Cp  is  the  fluid  specific  heat  at  constant  pressure 
(Cp  =  ayaCPa),  A  is  the  thermal  conductivity  of  the  fluid  and  ha  is 
the  sensible  enthalpy  of  the  species  a  in  the  fluid.  These  properties 
of  the  fluid  are  temperature  dependent.  The  JANAF  thermochemical 
tables  are  used  to  calculate  the  specific  heat  at  constant  pressure 
and  the  sensible  enthalpy  for  each  species  a  in  the  fluid: 

Cpa  =  TTj  (^la  +  a2oJ  +  a3aT2  +  <?4a^ 3  +  ^5aT4)  (12) 

Wa 


K-f 

Aref=o 


CpadT 


=  3L  (ai«T  +^t2  +  ^t3  +  + 


(13) 


where  aia,  a2a.  a3a>  a4a,  05a  are  the  JANAF  constants  [17].  The 
semi-empirical  formula  of  Wilke  (1950)  is  used  to  estimate  the 
multicomponent-fluid  thermal-conductivity  [16]: 


A  = 


^a  Arv 


(14) 


where 


0a(3  — 


l+(Xa/Xp  )1/2(lVp/Wa 


s\/4 


(8  +  8  Wa/Wp) 


1/2 


(15) 
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The  Sutherland  model  accounts  for  the  temperature  dependence  of 
the  thermal  conductivity  of  each  species  a  in  the  fluid: 


•( 


\  ^  ?o Xa  +  S\a 


'  oXa 


T  +  S; 


Xa 


(16) 


where  the  Sutherland-law  parameters  (A0,  T0x,  S^)  are  tabulated 
for  the  most  common  gases  in  [16]. 

2.2.  Electrode  model 

In  the  electrodes,  the  same  physical  phenomena  as  in  the 
channels  are  studied.  However,  the  governing  equation  for  each 
phenomenon  differs  from  that  in  the  channel  due  to  the  porous 
nature  of  the  electrodes.  The  momentum  conservation-equation  is 
formulated  in  the  form  of  Darcy’s  Law: 

v=  -—Vp  (17) 

Ji 

where  v  represents  the  superficial  permeation  velocity,  B0  the 
porous-medium  permeability,  p  the  pressure,  and  fi  the  gas- 
mixture  viscosity.  Eq.  (14)  is  similarly  applied  to  estimate  the 
viscosity  of  the  multicomponent  fluid: 

Xat^a 


*=£ 


(18) 


where  /xa  represents  the  viscosity  of  each  of  the  species  a  in  the 
gas  mixture  and  is  given  by: 

l+(/WMp)1/2(  Wp/Wa)1/4]2 


4>ai 8  = 


(8  +  8  Wa/Wp)1/2 


(19) 


The  viscosity  of  each  species  a  is  estimated  from  the  Sutherland 
law,  as  follows: 

3/2  „ 


P'a  —  P'oa 


iO|xa 


+  Si 


|xa 


[0|jLa  J  T  +  S|xa 


(20) 


The  Sutherland-law  parameters  /x0,  T0pu,  S ^  are  tabulated  in  [16]. 
Pressure  is  calculated  from  the  mixture  composition  using  the 
equation: 


P  =  2_J>a 


(21) 


where  pa  is  the  partial  pressure  of  species  a. 

The  equation  for  the  conservation  of  the  chemical  species  a  in 
the  gas  mixture  is: 


V  •  Na  =  0 


(22) 


where  Na  is  the  total  molar  flux  of  species  a,  given  by  the  Dusty 
Gas  Model  [12]: 


Vpa  _  1  ~  PaNp  Na  Da  B0 

—RT^  ~p  Di?F  +  +  RTD^Ill 

P  a  a3  al<  aK^ 


Vp 


(23) 


where  R  is  the  ideal  gas  constant,  T  is  the  fluid  temperature,  is 
the  effective  binary  diffusion  coefficient  of  species  a  in  species  (3: 

(24) 


0$  = 


a  3 


where  £  is  the  electrode  porosity,  r  is  the  electrode  tortuosity  fac¬ 
tor  and  Dap  is  given  by  Eq.  (6);  and  is  the  effective  Knudsen 
diffusion  coefficient  of  species  a  [18]: 


eff  _  £  dp  /  8 RT 
U aK  “  r  3  V  nWa 


(25) 


After  some  manipulation  of  Eq.  (23),  Na  may  be  expressed  as: 


=  -CJVP<*  +  Pet  +  V%*Pa 

(26) 

where: 

(27) 

RT(EP^VDfp  +  VDfK) 

p* 

=  —%rV 

Deff 

^aK 

(28) 

(29) 

3  ^  a 

Further  details  of  the  above  manipulation,  for  the  isothermal  case, 
are  given  in  [7].  The  physical  meaning  and  relative  relevance  of  r*, 
u£*  and  are  investigated  in  [8],  also  for  the  isothermal  case. 
Inserting  Eq.  (26)  into  Eq.  (22),  the  conservation-equation  for  the 
species  a  may  be  expressed  as  follows: 

-V  •  (r*vpa)  +  V  •  (t)P*Pa)  +  V  •  (C£*vpa)  =  0  (30) 

Reaction  source  terms  are  not  present  on  the  right  hand  side  of  Eq. 
(30)  because  the  electrochemical  reaction  is  modelled  as  a  super¬ 
ficial  reaction.  The  species  sources  or  sinks  in  the  electrodes  are 
given  by  Faraday’s  law  and  are  treated  as  boundary  conditions  at 
the  electrode-electrolyte  interface  (reaction  wall,  rw): 

Na,rw  =  (31) 

where  I  is  the  current  density  at  the  reaction  wall  ( Ja  or  Ic  depending 
on  the  electrode;  see  Section  2.3),  n  represents  the  number  of  elec¬ 
trons  to  convert  a  single  molecule  of  species  a,  F  is  the  Faraday’s 
constant  and  nrw  is  the  reaction-wall-normal  unit  vector.  Finally, 
principle  of  energy  conservation  is  applied  to  calculate  the  tem¬ 
perature  profile  within  the  electrodes.  The  following  assumptions 
apply  to  the  SOFC  electrodes:  (i)  local-thermal  equilibrium  may  be 
assumed  between  the  porous  matrix  and  the  fluid  flowing  through 
the  void  space  of  the  electrode,  as  reported  in  [23  ] ;  (ii)  the  volumet¬ 
ric  heat  sources  within  the  electrodes  are  neglected  since  they  are 
opaque  bodies  [20,24],  and  Joule  heating  is  only  considered  within 
the  electrolyte;  (iii)  dh  =  CpdT ;  and  (iv): 


V- 


(phv)  +  v  •  [  )  =  v  •  ( 


(32) 


From  Eq.  (9),  the  sensible  enthalpy  conservation  equation  in  the 
electrodes  turns  into: 


V- 


E1 


(WaNaCpa 


)T 


-  V  •  (A.effVT) 


=  uj  +  TV  • 


£ 


(WaNaCpa) 


~  2Jhc,V  '  (W«Na)] 


(33) 


where  Gpa  *  ha.  ,  Na  are  given  by  Eqs.  (12),  (13)  and  (26)  respectively; 
and  A,eff  is  the  effective  thermal  conductivity  of  the  fluid  and  porous 
medium  as  a  continuum,  modelled  according  to  the  upper  Wiener 
bound  [19]: 


Aeff=(l  -£)hs+£h 


(34) 


where  k  represents  the  thermal  conductivity  of  the  fluid  and  is 
given  by  Eq.  (14),  while  ks  is  the  thermal  conductivity  of  the  solid. 
The  heat  release  due  to  the  electrochemical  reaction  is  assumed 
to  take  place  on  the  anode  side.  Thus  in  Eq.  (33):  m  =  0  in  the 
cathode;  while  in  the  anode  the  molar  reaction-heat  release  is 
UJ'  =  -  242000  -  6.05(Tref  -  298)  [J  mol”1  ]. 
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The  effect  of  surface-to-surface  radiative  heat-exchange  is 
included  using  a  view-factor-based  radiation-model.  The  model 
calculates,  in  the  manner  described  in  the  Appendix,  the  heat  flux 
qrad  arriving  at  or  departing  the  surfaces  of  the  electrodes  from  or  to 
all  of  the  surfaces  in  the  domain;  this  flux  is  included  in  the  energy 
equation  for  the  electrodes  as  a  boundary  condition. 

2.3.  Electrolyte  model 


The  total  electromotive  force  of  the  above  electrochemical 
system  is  described  by  the  Nernst  equation.  However,  during  oper¬ 
ation,  the  real  voltage  of  a  cell  is  smaller  than  that  predicted  by 
Nernst  due  to  the  cell  resistances,  and  then  the  voltage  of  the  cell 
is  given  by 


V  £o  ,  ^ln/PH2V/(PQ2/pref) 

2F  \  Ph2o 


Pohm  -  (Pact  +  Pcon)a 


The  electrolyte  is  an  impervious  solid,  that  impedes  gas  flow. 
Nevertheless,  there  are  mass  and  charge  transport,  and  heat  release, 
due  to  the  anionic  (O2- )  current.  The  charge  transport  is  assumed  to 
be  one-directional  (normal  to  the  reaction  walls);  this  assumption 
is  valid  for  thin  electrolytes,  such  those  in  the  microtubular  anode- 
supported  cells  (^20  /im).  In  tubular  cells  the  reaction  area  in  the 
cathode  (Ac)  differs  from  that  in  the  anode  (Aa)  since  the  electrolyte 
has  a  finite  thickness.  Hence,  the  current  densities  at  both  reaction 
walls  (active  electrode-electrolyte  interfaces)  satisfy  the  following 
relationship,  which  ensures  charge  and  mass  conservation  in  the 
electrolyte: 

U  =  £/c  (35) 

Heat  transfer  through  the  impervious  electrolyte  is  only  possible 
by  means  of  conduction  and  radiation  mechanisms.  The  effect  of 
radiation  within  the  state-of-the-art  very  thin  electrolytes  has  been 
found  to  be  negligible  [20,24].  Conduction  is  thus  the  only  relevant 
heat-transfer  mechanism.  Thus,  from  Eq.  (9): 

-V  ■(kVT)  =  Q.  +  m  (36) 

where  zu  =  0  because  the  reaction  heat  is  considered  to  be  entirely 
released  in  the  anode;  and  the  volumetric  heat  sources  are,  in 
absence  of  radiation,  only  due  to  the  ohmic  heating  (Q.  =  Qohm),  also 
known  as  Joule  heating: 


-V  ■  (AV7)  =  Qohm 


(37) 


Vocv 

—  (Pact  +  Pc  on)c  (44) 

where  E°  is  the  standard  electrochemical  cell  voltage, 
E°  =  1.271  - 2.731  xlO_4I;  the  partial  pressures  of  the  involved 
species  are  referred  to  the  feeding  streams  compositions;  and  rjcon, 
Pohm*  Pact  stand  for  the  concentration,  ohmic  and  activation  losses. 
The  concentration  polarization  is  due  to  mass  transport  resistances 
in  the  electrodes;  it  is  calculated  as  follows: 


_  RT 

Peon,  a  —  2p 


In 


(  PH2PH20,rw 
yPH2oPH2,rw 


(45) 


_  RT 
Pc  on,c  —  2p 


In 


,rw 


(46) 


where  pa,rw  is  the  partial  pressure  of  the  species  a  at  the 
electrode-electrolyte  interface  (reaction  wall).  The  ohmic  polar¬ 
ization  is  mainly  caused  by  the  resistance  to  the  transport  of 
oxygen-anions  through  the  electrolyte;  it  may  be  estimated  as: 

8lej 

Pohm  —  (47) 

where  le  is  the  electrolyte  thickness.  Finally,  the  Butler-Volmer 
equation  provides  the  relationship  between  the  current  density  and 
the  activation  polarization  of  the  electrode.  The  most  general  form 
of  the  Butler-Volmer  equation  is: 


where: 

Qohm  =  %  (38) 

where  the  mean  current  density  in  the  electrolyte  is  estimated  as 
Ie  =  (/a  +  /c)/2;  and  the  electrolyte  anionic  conductivity,  cre,  is  given 
by: 

<7e  =  ae  exp  (-^)  (39) 

ae  and  be  are  constants  tabulated  in  [25]  for  the  common  electrolyte 
materials. 

The  effect  of  surface-to-surface  radiative  heat-exchange  is 
included,  as  for  the  electrodes,  as  a  boundary  condition  (see 
Appendix). 

2.4.  Electrochemical  model 

The  electrochemical  reaction  that  takes  place  in  a  H2-fed  SOFC 
consists  in  the  oxidation  of  the  hydrogen  in  the  anode  and  the 
reduction  of  oxygen  in  the  cathode,  releasing  a  water  molecule  as 


a  product  in  the  anode: 

H2^  2H+  +  2e-(anode)  (40) 

l/202  +  2e-  ->  02-(cathode)  (41) 

2H+  +  O2-  -*  H20(anode)  (42) 

The  overall  electrochemical  reaction  is  described  as  follows: 

H2  +  1/202^H20  (43) 


I  =  Io 


-  exp 


f  -oiFr] art 

^  RT 


(48) 


where  a  and  a  are  the  backward  and  forward  transfer  coefficients 
and  J0  is  the  exchange  current  density,  the  calculation  of  which  for 
both  electrodes  uses  the  following  experimental  correlations  [26]: 


fo,a  - 

=  Va 

/Ph2  \  /Ph2o> 
\Pref  /  V  Pref  7 

-0.5  , 

)  exp  l 

Fact,a  \ 

RT  ) 

(49) 

fo,c  z 

=  Yc( 

f  £act,c\ 

,  RT  ) 

(50) 

where  ya  and  yc  are  the  anodic  and  cathodic  pre-exponential  coef¬ 
ficients  and  Fact, a  and  Fact,c  are  the  anodic  and  cathodic  activation 
energies.  Eq.  (49)  applied  to  both  electrodes  provides  the  corre¬ 
sponding  activation  overpotentials:  pact,a  and  ijac t,c- 


3.  Numerical  details 


The  domain  considered  for  the  simulation  of  an  anode- 
supported  micro-tubular  SOFC  consists  of  five  concentric  and 
adjacent  cylindrical  sub-domains,  viz.  from  the  axis  outwardly: 
(i)  the  fuel  channel;  (ii)  the  anode;  (iii)  the  electrolyte;  (iv)  the 
cathode;  and  (v)  the  air  channel.  The  axial  symmetry  of  the 
tubular  geometry  simplifies  the  three-dimensional  problem  to  a 
two-dimensional  one.  Thus,  each  of  the  submodels  described  in 
Section  2  is  solved  in  its  corresponding  two-dimensional  mesh  and 
the  respective  solutions  are  coupled  as  required  through  boundary 
conditions.  Fig.  la  shows  a  schematic  of  the  domain,  representing 
an  anode-supported  micro-tubular  cell  inside  a  cylindrical  furnace. 
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Fig.  1.  Domain  diagram. 


The  anode  and  the  electrolyte  are  longer  than  the  oven  to  allow 
the  use  of  low- temperature  seals;  the  cathode,  i.e.  the  cell  active 
area,  is  centered  in  the  furnace  (Fig.  la.l)  for  the  cell  to  operate  at 
conditions  as  close  as  possible  to  isothermal  ones. 

The  channel,  electrode  and  electrolyte  models,  described  in  Sec¬ 
tion  2,  are  solved  numerically  using  the  finite-volume  method. 
Fig.  2  is  a  schematic  representation  of  the  solution  algorithm  (black- 
solid  lines)  and  the  coupling  paths  between  fields  of  adjacent 
meshes  (dashed  lines);  the  sub-indexes  stand  for  the  channel  inlet 
(in),  the  channel  outlet  (out),  the  heated  furnace  wall  (oven),  the 
electrode-electrolyte  reaction  walls  (rw),  the  channel  interface 
coupled  to  the  electrode  (Cl),  and  the  electrode  interface  coupled 
to  the  channel  (El).  In  the  channels,  the  continuity  and  momentum 
equations  are  solved  using  SIMPLE  [27].  In  the  continuity  equation 
(Eq.  ( 1 )),  the  volumetric  mass  source  term  (S)  is  zero  throughout  the 
domain.  However,  in  the  channel-electrode  coupled  boundary,  a 
correction  is  applied  to  ensure  continuity  through  both  domains.  In 
the  electrode  domain,  continuity  is  fulfilled  by  the  total  molar  fluxes 
given  by  the  DGM  (Na)  but  not  by  the  velocity  given  by  Darcy’s 


law;  continuity  through  the  coupled  boundary  is  thus  achieved  as 
follows:  (i)  va  =  x>Ei  and  (ii)  Sa  =  ^(N«, EiWa)  -  vapa.  Similarly, 

a 

in  Eq.  (3)  Sa  Ci  =  NaiEiWa  -  ^aPciya.ci-  The  noteworthy  feature  of 
this  algorithm  is  the  absence  of  coupling  between  the  tempera¬ 
ture  fields,  due  to  the  use  of  a  coupled-matrix  algorithm  to  solve 
simultaneously  Eqs.  (11),  (33)  and  (37).  Further  details  of  the  mass-, 
momentum-  and  species-fields  coupling  are  given  in  [7]. 

The  electrochemical  algorithm  is  solved  iteratively  and  this  is 
also  shown  in  Fig.  2  (dotted  lines).  The  cell  voltage,  V,  is  the  model 
input  and  is  assumed  constant  through  the  cell  active  area.  The 
mathematical  model  is  initially  solved  with  a  guessed  current  den¬ 
sity  field  at  the  cathode  reaction  wall;  the  current  density  for  the 
next  iteration  step,  Jc,new,  is  then  calculated  as  follows: 

*7ohm  +  (*7act  +  ^7con)a  +  (?7act  +  T]c  on)c 

fa, new  is  then  calculated  from  Eq.  (34).  The  new  estimated  current 
density  is  not  necessarily  constant  along  the  electrodes  since  the 


Fig.  2.  Numerical  algorithm.  Coupling  paths  are  shown  in  dashed  lines  and  the  electrochemical  loop  in  dotted  lines. 
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Table  1 

Cell  geometry  and  operating  conditions. 


Cell  geometry  [5] 

Fuel  channel  internal  radius 

1.2  mm 

Anode  thickness 

400  fxm 

Electrolyte  thickness 

20  fxm 

Cathode  thickness 

50  |xm 

Air  channel  external  radius3 

1  cm 

Fuel  channel  length 

10cm 

Anode  length 

10cm 

Electrolyte  length 

10cm 

Cathode  length 

1  cm 

Air  channel  length3 

7.5  cm 

Microstructural  parameters  [28] 

Anode  porosity 

41.4% 

Anode  tortuosity  factor 

3 

Anode  mean  pore  diameter 

1.2  |jim 

Anode  permeability 

7.65e-15  m2 

Cathode  porosityb 

50% 

Cathode  tortuosity  factor13 

3 

Cathode  mean  pore  diameter13 

1  fxm 

Cathode  permeability5 

7.65e-15  m2 

Operating  conditions  [5] 

Oven  set  temperature,  T0Ven 

1123  K 

Fuel  temperature  at  inlet,  Tin  fuei 

293  K 

Operating  pressure,  pout 

100000  Pa 

Fuel  composition,  xain  -►  H2:H20:Ar 

4.85:3:  92.15% 

Air  composition,  xajin  -*  02:N2 

21: 79% 

Fuel  inlet  (@room),  forced  convection,  Qin.fuei 

0.21  min"1 

Air  inlet,  natural  convection,  vin 

0. Oli  ms-1 

3  Personal  Communication. 

5  Unmeasured,  typical  value  assumed. 

cell  overpotentials  depend  on  the  species  and  temperature  distri- 

butions. 

4.  Results  and  discussion 

4.1.  Model  validation 

The  performance  of  the  mathematical  model  and  the  numeri¬ 
cal  algorithm  presented  in  Sections  2  and  3  has  been  evaluated  by 
simulating  the  behaviour  of  a  real  hydrogen-fed  anode-supported 
micro-tubular  SOFC  and  comparing  the  numerical  results  with  the 
corresponding  experimental  data.  The  experimental  I-V,  I-P  curves 
reported  by  Campana  et  al.  (sample  2)  in  [5]  are  the  selected  data 
to  carry  out  the  validation.  The  cell  geometry,  the  operating  condi¬ 
tions  and  the  electrochemical  model  parameters  are  summarized  in 
Tables  1  and  2.  In  the  absence  of  the  experimental-characterization, 
typical  values  for  the  cathode  microstructure  have  been  assumed. 

Table  2 

Electrochemical  and  thermophysical  properties. 

Electrochemical  parameters 

aa,  «a,  Uc,  «ca 

0.5 

Fact, a.  Fact,c  [26] 

120  kj  mol-1 

Yab 

le  + 11  Am-2 

Ycb 

le  +  9  Am-2 

Thermal  properties  [25] 

Anode  thermal  conductivity 

3  Wm1  K"1 

Electrolyte  thermal  conductivity 

2  Wm1  K"1 

Cathode  thermal  conductivity 

4  Witt1  K"1 

Electrolyte  emissivity 

0.4 

Cathode  emissivity 

0.4 

Furnace  emissivity 

0.8 

Electrical  properties  [25] 
ae 

85,000  Sm-1 

be 

ll.OOOK 

^exp  n  VSjm  ■  PeXp  O  Psim  • 

Fig.  3.  Experimental  data  (open  symbols)  [5]  vs  simulation  results  (solid  symbols) 
for  the  cell  in  Tables  1  and  2.  The  left  axis  is  the  cell  voltage  (V)  and  the  right  axis  is 
the  cell  power  ( P=IV ). 

This  assumption  will  not  bear  a  significant  influence  on  the  results 
since  the  concentration  overpotential  in  such  thin  cathode  is  neg¬ 
ligible  [5].  Numerical  results  are  plotted  in  Fig.  3;  where  the  x-axis 
represents  the  mean  current  density  over  the  cathode-active  area 
(/  =  A~ 1  fA  IcdAc );  good  agreement  with  the  experimental  curves 
is  shown. 

The  contribution  of  each  overpotential  (mean  value,  i)*  = 
A~ 1  J  rj*dA*)  to  the  total  voltage  loss  of  the  cell  given  in  Table  1  is 
shown  in  Fig.  4.  Again,  the  numerical  results  are  in  agreement  with 
the  experimental  ones  reported  in  [  5  ] ;  the  cathode-activation  resis¬ 
tance  and  the  anode-concentration  resistance  are  the  two  main 
causes  for  the  voltage-drop,  while  the  cathode-concentration  resis¬ 
tance  may  be  neglected. 

4.2.  Discussion 

As  indicated  in  the  Introduction,  SOFC  technology  must  still 
overcome  some  issues  before  becoming  widely  adopted  commer- 


Tlact,c  ®  ^acta  ®  hcon.c  A  ^con.a  ^  ^ohm 


Unmeasured,  typical  value  assumed. 
Fitting  parameters. 


Fig.  4.  Mean  overpotentials  at  different  mean-current  densities  for  the  cell  given  in 
Tables  1  and  2. 
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Fig.  5.  I-V-P  curves  for  the  cell  given  in  Tables  1  and  2  with  preheated  fuel  at 
different  fuel  flow-rates. 


Fig.  7.  Maximum  cell  temperature  at  different  mean-current  densities  for  the  cell 
running  with  preheated  fuel  at  different  fuel  flow-rates. 


dally.  Decreasing  the  operating  temperature  has  been  reported  as 
one  of  the  key  factors  to  face  those  issues,  since  lower  operating 
temperature  brings  about:  (i)  more  economical  cell  sealants,  thus 
reducing  the  overall  SOFC  costs;  (ii)  smaller  cell  thermal  stresses, 
thus  improving  the  cell  reliability;  (iii)  faster  start-ups;  and  (iv) 
smaller  diffusion  rates  of  the  electrodes  transition  metals  into  the 
electrolyte,  which  reduces  the  degradation  rate  of  the  cell  mate¬ 
rials  and  improves  their  durability.  Heat  management  thus  plays 
an  important  role  to  minimize  the  temperature  gradients  and  the 
hot  spots  within  the  SOFC.  The  temperature  distribution  is  not  easy 
to  determine  experimentally.  Therefore,  a  numerical  tool  becomes 
essential  to  optimize  the  heat  management  of  the  SOFCs.  Several 
numerical  works  on  the  temperature  distribution  in  SOFCs  may 
accordingly  be  found  in  the  literature.  For  instance,  Oulmi  et  al. 
[30]  present  a  numerical  study  of  the  temperature  distribution  in 
a  planar  SOFC  as  a  function  of  the  SOFC  configuration  (anode  or 
electrolyte  supported);  the  model  by  Jeon  et  al.  [31]  is  applied  to 


study  the  temperature  field  in  a  SOFC  and  its  components  along 
and  across  the  cell  at  different  loads;  Qu  et  al.  [32]  predict  the  tem¬ 
perature  distribution  within  an  anode-supported  planar  SOFC  with 
corrugated  bipolar  plates;  and  Serincan  et  al.  [29]  use  their  SOFC 
model  to  study  the  influence  of  the  preset  operating  temperature  on 
the  cell  performance  and  the  temperature  distribution.  In  all  these 
numerical  works,  the  authors  have  assumed  the  inlet  temperature 
of  the  feeding  gases  (air  &  fuel)  to  be  equal  to  the  preset  furnace- 
temperature.  This  assumption  is  often  needed  either  because  the 
experimental  test  rig  is  not  described  in  detail,  or  to  reduce  the  cal¬ 
culation  time.  However,  the  above  assumption  is  only  applicable  if 
the  gases  are  sufficiently  preheated  before  reaching  the  cell  active 
area;  otherwise  the  gases  may  indeed  cool  down  the  cell,  resulting 
in  counter-intuitive  cell  performance,  as  it  is  illustrated  below. 

Suzuki  et  al.  [4]  reported  a  significant  improvement  in  the  per¬ 
formance  of  their  anode-supported  microtubular  SOFC  as  they 
increased  the  fuel  flow.  This  effect  was  numerically  studied  by  Ser¬ 
incan  et  al.  [29],  who  attributed  the  performance  enhancement  to 
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Fig.  6.  Mean  overpotentials  at  different  mean  current  densities  for  the  cell  given  in  Tables  1  and  2  with  preheated  fuel  at  different  fuel  flow-rates. 
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Fig.  8.  Temperature  contours  (in  Kelvin)  for  the  cell  given  in  Tables  1  and  2  with  fuel  preheating  (Tin  fuei  =  1123 1<)  at  0.7  V  and  different  fuel  flow-rates.  The  cell  active  area  is 
indicated  with  black  straight  line. 


the  decrease  of  the  anode  concentration-overpotential,  remarking 
that  the  high  fuel  flow-rates  increase  the  power  output.  In  this 
section,  the  effect  on  this  outcome  of  the  widely  used  assump¬ 
tion  that  the  fuel  inlet  temperature  equals  the  oven  temperature, 
Tin, fuel  - ’Toven,  is  studied  for  the  cell  given  in  Table  1,  for  which 
the  fuel  is  usually  fed  at  room  temperature  to  the  rig  described 
in  Fig.  1.  The  cell  is  thus  simulated  operating  at  four  fuel  flow-rates 
(Qin, fuel  =  0-1  >  0.2,  0.3,  0  .41  min-1,  referred  to  room  temperature 
and  pressure)  and  at  two  fuel  inlet-temperatures,  viz  preheated 
flow  (Tjn  fuel =  ?oven  =  1123 1<)  and  room  temperature  (Tin  fuel  =  293  K, 

TOVen  =  li23K). 

Fig.  5  shows  the  I-V-P  for  the  given  fuel  inlet-flows  assuming 
preheated  fuel.  The  results  show  that  high  fuel  flows  improve  the 
cell  performance,  which  is  in  agreement  with  the  experimental  and 
numerical  results  of  Suzuki  et  al.  [4]  and  Serincan  et  al.  [29].  The 


results  plotted  in  Fig.  6  also  agree  with  Serincan’s  findings  [29],  i.e. 
the  improvement  of  the  cell  performance  when  increasing  the  fuel 
flow  is  mainly  due  to  a  decrease  in  the  anode  concentration  overpo¬ 
tential.  The  fuel  flow-rate  when  the  fuel  is  preheated  does  not  alter 
significantly  to  the  maximum  temperature  in  the  cell,  as  shown  in 
Fig.  7,  or  the  temperature  distribution,  as  shown  in  Fig.  8,  and  so  the 
temperature  dependent  overpotentials  ( rjact ,  ^ohm)  remain  almost 
constant,  as  shown  in  Fig.  6.  However,  Fig.  9  shows  the  I-V-P  curves 
for  the  cell  reported  in  [5]  at  different  fuel  flow-rates.  The  fuel  is  fed 
at  room  temperature  as  it  is  done  in  [  5  ] .  From  Fig.  9  it  is  clear  that  the 
cell-performance  improves  as  the  fuel  flow-rate  is  increased  from 
0.1  lmin-1  to  0.21  min-1.  However,  if  the  fuel  flow-rate  is  further 
increased  the  performance  of  the  cell  either  does  not  improve  (from 
0.2 1  min-1  to  0.3 1  min-1 )  or  it  even  deteriorates  (from  0.3 1  min-1  to 
0.4 1  min-1 ).  As  it  is  shown  in  Fig.  1 0,  large  fuel  flow-rates  decrease 
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Fig.  9.  I-V-P  curves  for  the  cell  given  in  Tables  1  and  2  at  different  fuel  flow-rates. 
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Fig.  10.  Mean  overpotentials  at  different  mean-current  densities  for  the  cell  given 
in  Tables  1  and  2  and  for  different  fuel  flow-rates. 


M.  Garda-Camprubi  et  al  /  Journal  of  Power  Sources  196(2011)  7290-7301 


7299 


■W.  =  293K 


* 

1 


0.8 ms-1.  It  should  be  mentioned  that  the  results  by  Suzuki  et  al. 
[4]  at  823  K  show  however  monotonic  improvements  in  the  cell 
performance  with  increasing  the  fuel  flow-rate.  This  however  does 
not  contradict  our  findings  regarding  this  cooling  effect,  since  this 
effect  becomes  less  relevant,  or  even  negligible,  at  lower  operating 
temperatures. 

Thus  the  positive  effect  of  the  increased  fuel  flow  on  the  con¬ 
centration  overpotential  may  be  countered,  depending  on  the  fuel 
feeding-temperature,  by  the  cooling  effect  of  the  flow.  Among 
other  corollaries,  this  highlights  the  importance  of  the  precise 
knowledge  of  the  operating/boundary  conditions  when  respec¬ 
tively  conducting  and  interpreting  experiments/simulations.  As  an 
example,  if  Tin  fuei  =  Toven  is  improperly  assumed,  then  the  effects 
of  the  induced,  and  neglected,  fuel-side  cooling  on  the  cell  per¬ 
formance  would  be  wrongly  attributed  to  the  electrochemistry 
parameters,  which  are  often  fitted  from  the  results. 


r/Acm'2 
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5.  Conclusions 


Fig.  11.  Maximum  cell  temperature  at  different  mean-current  densities  for  the  cell 
given  in  Tables  1  and  2  and  for  different  fuel  flow-rates. 

the  anode  concentration-overpotential  (as  expected),  but  the  acti¬ 
vation  overpotential  is  significantly  increased  due  to  a  cooling  effect 
shown  in  Fig.  1 1 ;  this  effect  is  often  ignored  in  numerical  simula¬ 
tions  where  the  assumption  Tin,fuei =  Toven  is  made.  Fig.  12  shows  the 
temperature  distribution  inside  the  test  rig  when  the  cell  is  oper¬ 
ating  at  0.7  V  without  external  fuel  preheating  for  the  three  fuel 
inlet- velocities  considered.  As  seen  in  Fig.  12,  the  fuel  is  heated  as 
it  flows  inside  the  furnace  from  the  inlet  to  the  active  area  of  the 
cell;  at  low  fuel  flow  rates  the  fuel  reaches  the  active-area  at  about 
the  preset  temperature  as  desired.  However,  at  higher  fuel-flow 
rates  the  path  from  the  inlet  to  the  active-area  is  not  long  enough 
for  a  full  preheating  of  the  gas.  Fig.  1 1  also  shows  that  the  cell  heat 
release  even  at  high  current  densities  may  not  be  enough  to  heat 
the  large  non-preheated  fuel  flow.  These  results  are  also  in  agree¬ 
ment  with  the  experimental  curves  reported  by  Suzuki  et  al.  [4], 
who  found  a  significant  improvement  of  the  cell  performance  when 
increasing  the  fuel  inlet-velocity  from  0.2  m  s-1  to  0.4  m  s-1  (873  K), 
but  no  major  changes  when  further  increasing  it  from  0.6  m  s-1  to 


This  paper  has  presented  a  comprehensive  CFD  model  of  the 
main  mass  and  heat  transfer  processes  taking  place  in  a  SOFC. 
The  model  uses  five  distinct  but  coupled  subdomains  for  the  two 
channels,  the  two  electrodes  and  the  electrolyte.  An  algorithm  for 
coupling  the  solutions  among  these  domains  has  been  developed 
and  presented. 

The  model  has  been  validated  by  comparison  with  experimental 
results  from  a  laboratory  anode-supported  tubular  cell.  Further,  the 
model  has  been  exploited  to  investigate  the  influence  of  the  mass 
flow  rates  and  thermal  fields  in  the  cell  performance.  It  has  been 
shown  that,  due  to  convective  cooling,  the  cell  performance  may 
deteriorate  for  larger  mass  flow  rates  if  the  feeding  temperature  is 
low. 


Acknowledgements 

This  work  is  supported  by  the  Science  and  Innovation  Ministry 
of  the  Spanish  Government  ( Ministerio  de  Ciencia  e  Innovation,  Gob- 
ierno  de  Espana )  under  project  ENE2008-06683-C03-03/CON.  The 


1123K 
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Table  Al 

Summary  of  surface-to-surface  radiation  in  the  air  channel. 


A- 

Electrolyte 

Cathode 

Furnace  wall 

Furnace  front-end 

Furnace  back-end 

A 

A 

A 

A 

A 

Electrolyte 

X 

X 

V 

V 

V 

Cathode 

X 

X 

V 

V 

V 

Furnace  wall 

V 

V 

V 

V 

V 

Furnace  front-end 

V 

V 

V 

X 

V 

Furnace  back-end 

V 

V 

V 

V 

X 
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Appendix  A.  Appendix 


For  calculating  the  radiative  exchange  among  surfaces,  these  are 
discretized  into  a  number  of  elementary  surface  elements.  In  this 
work,  these  elements  are  also  the  relevant  faces  of  the  computa¬ 
tional  mesh  used  for  the  finite-volume  method.  The  radiative  flux 
departing  a  surface  element  i,  qra d(I-,  is  given  by  [22]: 


As  indicated  above,  surface-to-surface  radiation  in  the  cell  chan¬ 
nels  is  the  only  radiation  mechanism  relevant  in  the  operation  of  a 
solid-oxide  cell.  In  the  present  model,  the  appropriate  heat  fluxes 
are  calculated  using  a  view-factor  method,  and  inserted  in  the  heat- 
transfer  equations  for  the  electrolyte  and  electrodes  as  boundary 
conditions  at  the  electrode-channel  and  electrolyte-channel  inter¬ 
faces.  This  Appendix  provides  some  additional  details  on  the  model 
implemented. 

In  the  geometry  analysed  in  this  paper  [5],  the  fuel  channel 
consists  of  a  long  thin  cylinder  (2.4  mm  diameter  and  100  mm  in 
length),  made  of  a  porous  material  (anode).  The  radiative  heat  flux 
at  a  given  point  of  the  fuel  channel  can  thus  arise  only  from  irra¬ 
diation  from  other  locations  at  the  anode  surface.  Considering:  (i) 
axial  symmetry;  and  (ii)  view  factors  vanish  quickly  in  the  axial 
direction  for  such  a  thin  cylinder;  the  surface-to-surface  radiation 
in  the  fuel  channel  is  neglected. 

The  air  channel  consists  of  the  void  space  between  the  cylin¬ 
drical  furnace  and  the  concentric  cell.  The  cell  external  radius  is 
1.62  mm,  in  the  non-reacting  zone,  and  1.67  mm  where  the  cath¬ 
ode  is  printed  (see  Fig.  1).  The  furnace  is  10  mm  in  radius.  Both 
ends  of  the  air  channel  are  filled  with  a  thermal-insulating  porous- 
material,  retaining  the  heat  in  the  furnace  while  allowing  the  air 
go  through.  Surface-to-surface  radiation  thus  takes  place  between 
several  pairs  of  surfaces,  as  summarized  in  Table  Al  and  depicted 
in  Fig.  Al. 

From  the  above  discussion  it  is  concluded,  that  there  are 
five  surfaces  involved  in  the  surface-to-surface  thermal-radiation 
exchange:  (i)  the  electrolyte  external  surface;  (ii)  the  cathode  exter¬ 
nal  surface;  (iii)  the  furnace  heated  wall;  (iv)  the  furnace  front-end; 
and  (v)  the  furnace  back-end.  These  are  all  included  in  the  radiation 
model. 


9rad,i  —  ^rad,i 


Fi-j  Qr  ad,  j 


(Al) 


where  the  subindex  j  stands  for  each  of  the  other  surface  elements 
that  are  viewed  from  element  i;  £rad  is  the  emissivity  of  the  surface, 
Eb  is  the  blackbody  emissive  power,  E0  is  the  incident  radiation 
entering  or  leaving  the  enclosure  through  an  opening,  and  Fz_j  is 
the  view  factor  between  each  two  surface  elements,  At  and  Aj.  The 
view  factor  is  defined  and  calculated  as: 


FH  = 


cos(0,)cos(6U  -(n, ••?)(%•  r) 


7rrz 


-Aj  = 


jrr* 


(A2) 


where  6[  and  0j  are  the  angles  between  the  surface  normal  vectors 
and  the  line  connecting  A*  and  Aj  (of  length  r),  as  shown  in  Fig.  Al. 
This  computation  of  view  factors  is  performed  only  once,  at  the 
beginning  of  the  calculation. 

The  radiative  heat  flux  at  each  surface  element  is  calculated  as 
the  solution  of  the  system  of  equations  given  in  Eq.  (Al ) 

Regarding  the  emissivities  of  the  different  surfaces,  the  cell  elec¬ 
trolyte  is  made  of  Yttria  Stabilized  Zirconia  (YSZ),  the  emissivity  of 
which  is  uncertain.  Ferreire  et  al.  [33]  measured  the  YSZ  emissivity 
at  1273  K,  finding  a  value  of  0.25;  Alaruri  et  al.  [34]  calculated  an 
effective  emissivity  of  0.4  for  the  YSZ  samples  without  aluminium 
impurities  at  1123 1<;  Sully  et  al.  [35]  measured  the  emissivity  of 
the  pure  zirconia  at  1073K,  obtaining  a  value  around  0.32.  The 
electrolyte  emissivity  is  set  to  0.4  in  this  work. 


Furnace  wall 


Fig.  Al.  Computational  domain  for  view-factor  calculation  and  surface-to-surface  radiation  model. 
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The  cell  cathode  is  made  of  a  composite  material,  consisting  of 
Lanthanum  Strontium  Manganite  and  Yttrium  Stabilized  Zirconia 
(LSM-YSZ).  For  a  similar  material,  an  emissivity  of  0.4  has  been 
reported  [36],  and  is  used  in  this  model. 

The  material  of  the  furnace  enclosure  is  not  known,  and  its  emis¬ 
sivity  is  thus  assumed  to  be  0.8. 
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